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We generalize the Kuramoto model for coupled phase oscillators by allowing the frequencies to 
drift in time according to Ornstein-Uhlenbeck dynamics. Such drifting frequencies were recently 
measured in cellular populations of circadian oscillator and inspired our work. Linear stability 
analysis of the Fokker-Planck equation for an infinite population is amenable to exact solution and 
we show that the incoherent state is unstable passed a critical coupling strength K c (pj, 07), where 7 is 
the inverse characteristic drifting time and <7f the asymptotic frequency dispersion. Expectedly K c 
agrees with the noisy Kuramoto model in the large 7 (Schmolukowski) limit but increases slower as 7 
decreases. Asymptotic expansion of the solution for 7^0 shows that the noiseless Kuramoto model 
with Gaussian frequency distribution is recovered in that limit. Thus varying a single parameter 
allows to interpolate smoothly between two regimes: one dominated by the frequency dispersion 
and the other by phase diffusion. 

PACS numbers: 05.40.-a, 05.45.Xt, 89.75.-k, 87.19.Jj 
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I, INTRODUCTION 

Synchronization phenomena occur recurrently in phys- 
ical, chemical and biological systems. Few noteworthy 
examples include superconducting currents in Josephson 
junction arrays pill], emerging coherence in populations 
of chemical oscillators 0], or the accuracy of central cir- 
cadian pacemakers in insects and vertebrates 0, 0, @. 
The latter serve as biomolecular time-keeping devices, 
which most organisms have evolved to coordinate their 
physiology and metabolic activities with the geophysical 
light-dark and temperature cycles @]. 

The present work was motivated by recent experiments 
in mammalian cell cultures in which the levels of pro- 
teins implicated in the circadian (with ~ 24hr periods) 
clockwork were monitored using fluorescent reporters 
0, 0, El- It was demonstrated that individual cellu- 
lar oscillators generate self-sustained rhythms in protein 
abundance and that populations can be synchronized by 
treatment with a short serum shock or light pulse. Impor- 
tantly frequencies of individual oscillators are not strictly 
constant but drift in time (see for example Fig. S2 in [lfj| 
showing the results in the zebrafish). Several chemical 
kinetics models thought to capture the biochemistry re- 
sponsible for generating oscillations in living cells were 
shown to exhibit oscillatory instabilities and limit-cycles 
0,0]. Experimental and theoretical evidence therefore 
supports a description of oscillator populations in terms 
of phase variables. 



Our understanding of the onset of collective syn- 
chronization in coupled nonlinear oscillator models has 
greatly benefited from a large body of work on the Ku- 
ramoto model [H Q [H EH, Ej 
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describing the phase dynamics in a set of weakly coupled 
identical non-linear oscillators. Here, (fi(t) represent the 
phase of the z-th oscillator at time t, and are white 
noise sources with expectation and covariance 

E&(i)] - , 
Cav[fc(a),&(i)] - 2D8 ij 8{s-t) . 

The frequencies fi are static and taken from a distri- 
bution g(f) symmetric around fif . Eq. is an effec- 
tive model for the phase degrees of freedom in a popula- 
tion of limit-cycle oscillators and assumes a regime where 
phase and amplitude dynamics decouple. The parame- 
ter K measures the strength of the all-to-all coupling. 
Most exact results are given for the coupling function 
U(ifij — ipi) = sin(ipj — ipi). More general interactions 
lead to much gre ater analytical complexity and were in- 
vestigated in jlH. ITfll |. Critical properties of the model 
are conveniently studied using the complex order param- 
eter R(t)e**® = ^E^=ie iw(t) so that collective syn- 
chronization occurs when Roo = lim.x->oo y Jq " R(t) dt 
remains positive in the infinite population limit. For 
the sine coupling model a bifurcation occurs at K c = 
2// D'^+f 1 9 if + M/)^/ a t which the incoherent desyn- 
chronized state i?oo = becomes unstable and a macro- 
scopic number of oscillators phase lock to the average 
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phase ip(t) = fift HIH Il- 
eal Kuramoto result K c 



For D — > the classi- 
231 is recovered. Be- 
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low the critical coupling, the incoherent state is linearly 
stable when D > |22| but only neutrally stable when 
D = with i?(f) still decaying to zero [T^ |. 



II. THE MODEL 

To study the effects of the reported drifts on collective 
synchronization, we generalize the Kuramoto model by 
introducing a second time scale I/7 (besides I/07) char- 
acterizing the frequency drifts. The frequency dynamics 
is formulated as an Ornstein-Uhlenbeck (O-U) process 
while the phases are coupled following the canonical all- 
to-all sine interaction. The model for N oscillators reads 



time properties. For fixed 07, we anticipate that syn- 
chronization should be hardest for strictly static oscil- 
lators (7 = 0). This case corresponds to the noiseless 
D = Kuramoto model with Gaussian frequency dis- 
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so 



persion g[f) = J 

that K c — 2yj2/-K <j f . As the frequency dynamics loses 
stiffness (when 7 increases), we expect the synchroniza- 
tion threshold to be facilitated by the frequency drifts. 

We study the infinite population model N — * 00 by 
formulating a Fokker-Planck equation for the time de- 
pendent joint density p((p, /, t). The all-to-all interaction 
term 
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^awfeit) - <pi(t)) = KR{t)sin(ip(t)-(pi(t)) (4) 
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Ut) = - 7 (/i(t) -fi f ) + Vi (t) , 
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where [if is the average frequency chosen identical for 
each oscillator We assume that the rji are independent 
and identically distributed white noise sources with 

Efa(*)] = , 
Cov[r]i(s),r]j(t)] = T] 2 5ijS(s - t) . 

The solution for fa is a Gaussian process with mean and 
covariance 
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= ±Ls~(e-rt- 
7 t»i rf 



*l _ e -7(*+s) 
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In the following we use as independent parameters 
the asymptotic frequency dispersion o 2 = |- and the 
damping 7, which are in principle both accessible exper- 
imentally. Then, Cov[/i(s), fi(t)] — ► 2o 2 /j6(t — s) when 
ft 3> 1 . To remind the significance of this regime we 
note for K = the phases (pi(t) also follow Gaussian 

2er 2 

processes with Cov[ipi(t), >pi(t)] — > — asymptotically 
for jt 1. Because of the linear time dependence, this 
regime (Schmolukowski) describes phase diffusion with 

constant D = —. 

7 

Although it is a priori unclear whether this model 
exhibits a bifurcation, we expect that the large 7 be- 
havior reminiscent of phase diffusion will converge to 
the Kuramoto model (Eq. ^| with a frequency distri- 
bution given by <?(/) = S(f — fif) and white noise 
strength D = cry/7, and thus exhibit a bifurcation at 
K c = 2D. However the small 7 behavior is less ob- 
vious since we are simultaneously concerned with long 



has well known mean-field character and can be replaced 
for N — > 00 by 
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We obtain 
dp 
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d 2 p dp 
dp dip 
d(c(p, ip)p) 



f t = rff^js -f% + 7(/ " + IP (6) 
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This is the known expression for an O-U process 
augmented by a phase coupling involving c(j>, ip) = 
J n d9 fdgp(O,g,t)sm(0 — <p), which makes Eq. non- 
linear as a consequence of Eq. 03 



III. STABILITY ANALYSIS 

We next discuss the linear stability of the incoherent 
stationary solution po(f, /) = A/" M/j(T/ (f)/(2n) in first or- 
der. For reasons that will become clear, we factorize a 

1 ll 

term Jy^i a f(f) off the perturbation and write 

p(<P,f,t) =p (f)+Afy f 2 rTf (f)e( i p~ f i f t,aJ 1 (f~v f ), 1 t) , 

where e((p, /, t) is a small perturbation expressed in a 
rotating frame using rescaled frequency and time vari- 
ables. By plugging this ansatz into Eq. El we obtain the 
linearized problem || = Ce + 0(e 2 ) where C is the linear- 
operator 
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dgAf^ 2 (g)e(9,g,t)cos(6-ip) . 
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Decomposing s as a Fourier series in ip, e(<p, /, t) 
EnL-oo £ n(f, t)e~ mv , we obtain for the coefficients e n 



dt 



d 2 e n 

dp 



4 7 
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-*H«l^<i 2 (/) J dgN^ 2 (g)e n (g,t) 



(7) 



We notice that the first term representing the frequency 
dynamics resembles the harmonic oscillator plus a com- 
plex part, which can be removed by applying the trans- 
lation operator U$ defined by (Uqj\(x) = f(x — 9). We 

note that C n = U 2in ^£ n U_ 2in ^, where we have set 
a = (07/7) 2 and the operator C n — d 2 + \ — n 2 a — jf 2 
is self-adjoint on L 2 (R) and has pure point spectrum 
£(£„) = {X n( = -£-n 2 a : I = 0, 1, 2, . . .}. Its eigen- 
functions are given in terms of the Hermite functions j2J| 



H (x) 



and 



H t (x) = (2 e £l V^)~ 5 (-l)V 
as follows: 



= 1,2,. 



C n <S> f = \ nf <5>f 



Mf) = 2- 1 / 4 H e (2- 1 / 2 f) 



Therefore {U 2in ^§e '■ £ = 0, 1, 2, . . . } forms an orthonor- 
mal family which diagonalizes C n . Notice that the largest 
eigenvalue for each n is A„o = —n 2 a. 

Linear stability follows directly except for \n\ ^ 1 and 
K > 0. Indeed, for n = 0, we find A o = with cor- 

1 /2 

responding eigenfunction is $o = A^ j . However, this 
function lies outside the space of relevant perturbations 
because the normalization of p, J Q V dip J df p(ip, f,t) = 

1 /2 

1, requires orthogonality of Af ' x and eo(f,t) through 

J Nq^i Eo(f ,t)df = 0. Subsequent eigenvectors have 
negative eigenvalues. For all other \n\ ^ 1 the cou- 
pling term in Eq. vanishes and the incoherent state 
p(tp,f,t) — po(f) is linearly stable as a consequence of 
the strictly negative spectrum of C n . The same holds for 
all n in the absence of coupling K = 0. 

For the remaining case n — ±1 and K > 0, we notice 
that the coupling term in Eq. also vanishes for all di- 

1/2 

rections orthogonal to M \ , leaving a one-dimensional 
space that could develop an instability. We write the 
eigenvalue problem for Eq. implicitly as 

£„£„ + 5x|n| ^( £ «:-^0,l 2 )-^0,l 2 = ^ £ n ■ 

Using the resolvent equation 

(A — L n ) 1 = (A — U 2m ^£ n U _2m A /o) 1 

— u 2in ^(\ — C n ) 1 u_ 2in ^ 
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FIG. 1: Behavior of K c as a function of 7 as given by 
Eq. El with a — (07/7) 2 (continuous line). Dashed line rep- 
resents the Kuramoto model with identical frequencies and 
K c — 2D — 2<r 2 f/j. The 7^0 limit reproduces the 
7 = mode l wit h Gaussian frequency dispersion and gives 
K c /a f = 2^/2/^ = 1.5957 (dotted line). 



we obtain 



1/2 



K 

2l' 



) 2^ T — 

3=0 



"J 



where we have used the spectral decomposition C n f — 
'!-,,. / <!•,,. 

To find the critical coupling K c above which the inco- 
herent state becomes linearly unstable, we need to mon- 
itor when the largest eigenvalue crosses the imaginary 

1 /2 

After projecting onto M 01 , simplifying the fac- 
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tors (e n ,Af Q x ) on both sides of the equation, and setting 
A = we find an equation for K c : 



2_ _ f fa^U^g _ (-0)' 
K, ^ -x- 6 ^ 



3=0 



3=0 



e a a- a / t°- 1 e-*dt = e a a- a 7 (a,a) , (8) 

/Q 

where 7(0, x) is the lower incomplete T-function. 

The behavior of K c together with the Kuramoto model 
asymptotes for 7 — * 00 and 7 — > + limit are shown in 
Fig. It is noticeable that we find a bifurcation for all 
values of 7. K c strictly decreases from a finite 7 = 
as 7 increases, asymptotically behaving as K c — 2a 2 /■y. 
The analytical result thus supports the following picture: 
for small 7, the dominant source of fluctuations against 
which the coupling must work to achieve synchroniza- 
tion is the (Gaussian) frequency dispersion. As 7 in- 
creases while af is kept fixed, faster frequency drifts help 
synchrony by preventing individual oscillators with de- 
tuned frequency to stay out of tune for too long. Indeed 
with drifting frequencies every individual oscillators fluc- 
tuates around the mean frequency fXf with a time scale 
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7 _1 . In the large 7 regime, the effective frequency dis- 
persion vanishes and the coupling force needs to synchro- 
nize noisy but otherwise identical frequency oscillators. 
As predicted by the phase diffusion limit, the effective 
white noise strength D and hence K c decrease as 7 _1 . 

We now discuss the asymptotic regimes in detail: the 
small 7 limit follows from reverting to the original vari- 
ables and using the asymptotic expansion of 7(0, a) (us- 
ing Stirling's formula and 0: 6.5.3, 6.5.22, and 6.5.35). 
We obtain in the limit 7 — * + 

2 {^j) ~ V2 + 3^ + ^r^ + ° (7) - 

This proves that the model continuously interpolates to 
the noiseless (D = 0) model and that the 7^0 recovers 
the 7 = transition p redic ted in the original Kuramoto 
model at K c /af — 2^/2/tt. In the opposite regime 7 — > 
00 (thus a -> 0) we find (0 6.5.12, 13.1.2) 

a~°e°7(a,a) = a~ l M(l, 1 + a, a) - a~ x (l + G(a)) , 

where M(-, •, •) is the confluent hypergeometric function. 
This leads K c ~ 2aj/~f + C( 7 ~ 2 ) and hence proves the 
convergence to the white noise model (Eq. with D — 

Finally, we mention a generalization that includes a 
white noise source in the phase equation (as in Eq. in 
addition to the correlated frequency fluctuations. This 

Ft 2 

leads an additional diffusion term —D-q^i in Eq. Fol- 
lowing the steps above readily extends Eq.|H|to 

|^ = e a a- ( - a+b '>-f(a + b,a) , 
K c 

where b = D/j, with similar qualitative behavior. In 
particular, K c asymptotes to 2(0^/7 + D) for large 7 
and has finite 7 — > limit. 

IV. NUMERICAL SIMULATIONS 

We have performed numerical simulations of Eq. [3 to 
explore the behavior of R(t) (see Eq. EJl and in par- 
ticular Roc in function of the reduced coupling K. r = 
(K — K c )/K c . To verify the analytical results and study 
the scaling Roc = kK@ above the bifurcation, we simu- 
lated a finite number of oscillators using the exact so- 
lution for the frequency part, leading to t he updates 
fi(t +dt) = fi(t) e"^ 4 + M/(l - e" 7dt ) + 77 <r f VI - e-^ dt 
and ipi(t +dt) = ip^t) + (/»(*) + J£ £\ sm(ipj - ipi))dt 
where 77 is a Gaussian random number. We used Eq. 21 
to compute R(t) and transients were removed by waiting 
until the solutions from two different initial conditions 
(Pi(t = 0) = and ipi(t — 0) taken randomly converged 
to the same trajectory. The steady state value Roc was 
subsequently estimated by averaging R(t) over time. 

Fig. El fully supports the analytical solution and also 
indicates that the behavior of K c above the bifurcation 
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(K - K c )/K c 

FIG. 2: Numerical simulation of Eq. |5] Estimation of Roc 
was obtained using finite size scaling for systems of sizes N = 
5000, 10000 and 20000. Eq. El was used for K c to set the 
reduced coupling K r = (K — K c )/K c . Values for 7 were 4(o), 
3(A), 2.5(+), 2(x), 1.8(o) and l-6(v) and 07 = 1. In each 
simulation, 10 5 time steps of size dt — 0.01 were performed. 
We verified that the dependence in the step size was weak. 
Inset: 1/N finite size scaling for 7 = 1.8. K r = is the 
solid line, smaller (resp. larger) K r are below (resp. above) 
K r — 0. The extrapolated value for 1/N = is used in the 
main panel. 



depends only weakly on 7 over the simulated range. To 
inspect more closely whether Roc ~ \fK~r as in the Ku- 
ramoto model, we used refined spacing and larger sizes 
in the vicinity of K r = + . As shown in Fig.UH the sim- 
ulations are compatible with an exponent (3 = 0.5, the 
slightly higher exponents probably reflect a finite size ef- 
fect. On the other hand k correlates negatively with 7 
which is visible in both Figs [2] and 



V. DISCUSSION 

We have extended the Kuramoto model to frequencies 
which can drift in time following Ornstein-Uhlenbeck dy- 
namics. The net effect is that the white noise source in 
Eq. H is replaced by colored noise (with a Cauchy dis- 
tributed power spectrum), hereby adding a new time 
scale describing memory or frequency stiffness to the 
problem. Apart from mean field coupling among the 
phases which introduces a non-linearity, the stochastic 
phase and frequency dynamics lead to a linear Fokker- 
Planck operator which can be solved. Consequently the 
linear stability of the incoherent state can be addressed 
analytically using spectral calculus. The expression for 
the critical coupling above which macroscopic phase co- 
herence emerges can be resummed and expressed in terms 
of incomplete T-functions. Asymptotic expansion for 
small and large 7 shows that the full model continuously 
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0.01 0.02 0.05 0.10 ' 2.0 ' 3.0 ' 4.0 

(K - K c )/K c 7 

FIG. 3: Critical behavior above the bifurcation using step 
dK r = 0.01. Here, 7 is 4(o), 3(A), 2(+), 1.6(x), 1.2(o) and 
0.8(v). Notice the log-log scale to emphasize the power law. 
Lines are fits to Roc = kK^ Systems of sizes N = 10000, 
20000 and 50000 were simulated with the same parameters 
and same scaling procedure as in Fig. [2] Right panels show 
parameter estimates from the left panel. 



depending on the proximity of individual frequencies to 
the population mode. As 7 increases, (while crj is held 
fixed) the frequencies lose their stiffness which results in 
a reduction in the critical coupling K c needed for syn- 
chrony. Finally for very rapidly drifting oscillators (large 
7) cancel out the frequency distribution and generate an 
effective white noise source acting on the phases of oth- 
erwise equal frequency oscillators. At the same time the 
locked and incoherent subgroups become indistinguish- 
able. For intermediate 7, our numerical simulations in- 
dicate that the model belongs to the same (3 = 0.5 uni- 
versality class as the Kuramoto model. 

Because of analytical tractability and few parameters 
we expect this solution to be relevant for oscillatory sys- 
tems in the presence of complex noise sources. Such cases 
include populations of neural oscillators or biochemical 
oscillators where bioluminescence recordings have shown 
how intracellular noise sources generate frequency disper- 
sion through drifts. 



interpolates between two limits of the original Kuramoto 
model: one dominated by noise (large 7) and the other by 
the frequency dispersion (small 7). Therefore, the cou- 
pling force must counteract different sources of fluctua- 
tions to induce collective synchrony in drifting frequency 
oscillators, depending on the regime set by 7. Specifi- 
cally, for slowly drifting (small 7) frequencies the model 
approaches the noiseless model (Eq. with D = and 
g(f) = Ap,, CT/ (/)) where the coupling splits the popula- 
tion into distinct locked and incoherent sub-populations, 
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